1. Field of the Invention
The present invention relates to a back projection method in image reconstruction, and particularly to a method for arranging data allocation and computation sequence to perform calculation of back projection in image reconstruction fast.
2. Description of the Related Art
An internal structure of a human body can be obtained noninvasively in the form of tomogram by CT (Computed Tomography) apparatus. Recently, the internal structure can be obtained as a three-dimensional structure by using a plurality of tomograms with the progress in the apparatus.
X-ray image is a projected image, and CT image can be reconstructed from sectional images by back projecting the X-ray projection data photographed from a plurality of angles. This basic principle is described mathematically in “Uber die bestimmung von funktionen durch ihre intergralwerte langsgewisser mannigfaltigkeiten (on the determination of functions from their integrals along certain manifolds)”, J. Radon, Berichte Saechsische Akademie der Wissenschaften, vol. 29, pp. 262-277, 1917, and various calculation methods have been invented since then.
Currently, the main reconstruction method is back projection (hereinafter referred to as BP method) which is one of the methods of performing a back projection process for each two-dimensional section. The BP method is described in detail, for example, in “Principles of Computerized Tomographic Imaging (Classics in Applied Mathematics, 33)” by Avinash C. Kak, Malcolm Slaney: ISBN: 089871494X, pp. 60-75. There are developed BP methods such as filtered back projection (hereinafter referred to as FBP method) and convolution back projection-(hereinafter referred to as CBP method).
As an enormous amount of time is needed in the calculation of CT reconstruction, exclusive hardware is constructed and image processing hardware is used in the related art, in order to speed up the calculation. For example, such a technique is described in “OS Towards a Unified Framework for Rapid 3D Computed Tomography on Commodity GPUs”, Fang Xu and Klaus Mueller, Stony Brook University, Oct. 29, 2003, and “Hardware-Accelerated High-Quality Reconstruction of Volumetric Data on PC Graphics Hardware”, Markus Hadwiger, 2001.
A method for generating a CT image according to the related art is described below in brief. FIG. 32 is a flowchart showing the outline of a calculation procedure in the CT image generating method according to the related art. FIG. 33 is a schematic view for explaining a standard CT image generating method of a CT apparatus in the related art. FIG. 34 is an explanatory view showing a process for generating a projected image from an original X-ray image with the CT apparatus in the related art. FIG. 35 is an explanatory view showing a process of reconstructing the original image from the projected image.
X-ray image is an image projected from a single angle, and is a one-dimensional image. Accordingly, as shown in FIG. 33, two-dimensional image is synthesized by obtaining one-dimensional images from a plurality of angles ranging from 0° to 180°, while an X-ray generator is rotated around a subject.
As shown in FIG. 34, original image taken by X-ray photography is projected (mapped) from orthogonal coordinates (x,y) to rotating coordinates (r,θ), thereby generating a projected image (Sinogram). The projection is performed by Radon transformation represented by the expression 1, so that transformed projected image data is input to a computer (step S51 in FIG. 32).p(r,θ)=∫f(x, y)ds  [Expression 1]
Then, before performing the back projection, a filtering process is performed in order to prevent the image from being blurred by numerical error (step S52). This process will be referred to as a filtering portion or a convolution portion. Then, back projection is performed to the filtering-processed projected image (step S53).
As shown in FIG. 35, in the back projection, each projected image is additively projected onto a reconstructed image line-by-line in accordance with the angle from which the projected image is obtained. This process is referred to as a back projection portion. The back projection (inverse Radon transformation) is performed by the expression 2, and the convolution calculation is performed by the expression 3.
                              f          ⁢                                          ⁢                      (                          x              ,              y                        )                          =                              1                          2              ⁢              π                                ⁢                                    ∫              0              π                        ⁢                          q              ⁢                                                          ⁢                              (                                                                            x                      ⁢                                                                                          ⁢                      cos                      ⁢                                                                                          ⁢                      θ                                        +                                          y                      ⁢                                                                                          ⁢                      sin                      ⁢                                                                                          ⁢                      θ                                                        ,                  θ                                )                            ⁢                              ⅆ                θ                                                                        [                  Expression          ⁢                                          ⁢          2                ]                                {                                                                              q                  ⁢                                                                          ⁢                                      (                                                                  r                        ′                                            ,                      θ                                        )                                                  =                                  ∫                                      h                    ⁢                                                                                  ⁢                                          (                                                                        r                          ′                                                -                        r                                            )                                        ⁢                                                                                  ⁢                    p                    ⁢                                                                                  ⁢                                          (                                              r                        ,                        θ                                            )                                        ⁢                                          ⅆ                      r                                                                                                                                                                h                  ⁢                                                                          ⁢                                      (                    r                    )                                                  =                                                      1                                          2                      ⁢                      π                                                        ⁢                                      ∫                                                                                          ω                                                                    ⁢                                                                                          ⁢                      exp                      ⁢                                                                                          ⁢                                              (                                                  jω                          ⁢                                                                                                          ⁢                          r                                                )                                            ⁢                                              ⅆ                        θ                                                                                                                                                    [                  Expression          ⁢                                          ⁢          3                ]            
On the other hand, in the CT image generating method, real number calculation is frequently performed, and a large amount of time is required in the calculation for reconstructing a tomogram after taking CT sinogram images. Particularly, a large amount of time is required for real number interpolation calculation in back projection calculation. Therefore, some methods are proposed for increasing the speed of back projection in the reconstruction of a tomogram.
FIG. 36 is a flow chart for explaining an acceleration algorithm for back projection calculation disclosed in JP-A-2003-38480. In this calculation method, interpolated stretched data to which the back projection data is stretched for m times is calculated by interpolation before back projection. Then back projection is performed onto a tomogram reconstruction region. At the time of back projection, interpolated stretched data is back projected, being specified by multiplying the projected coordinates of each point which is to be reconstructed for m times.
In the back projection calculation, interpolated stretched data to which the back projection data is stretched for m times is calculated by interpolation before back projection. Then the process time for calculation of angle θ is shortened, which is in the innermost loop among a triple loop calculation concerning coordinate x, coordinate y and angle θ. That is, to obtain a projected value p for an angle θ, interpolation is not performed on the spot. An approximate value which is obtained beforehand by the interpolated stretched data is used so that the processing time for calculation (step S611) can be shortened (e.g. see JP-A-2003-38480).
In the back projection calculation in the CT image generating method according to the related art, calculation time is shortened by skipping the interpolation calculation of step S611 in FIG. 36. However, it is not sufficient for speeding up the calculation time, because real-number calculation for coordinate calculation in step S609 in FIG. 36 is performed in the innermost loop of the triple loop.